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Ion  thruster  plume  interaction  has  been  studied  extensively  in  recent  years.  While 
most  existing  plume  models  have  focused  on  charge-exchange  ion  interactions  with  the 
spacecraft  and/or  plume-induced  contamination,  few  studies  have  examined  the  detailed 
physics  near  the  thruster  exit.  In  particular,  the  ion  beam  neutralization  process  and  the 
characteristics  of  the  neutralizing  electrons  are  not  well  understood.  This  paper  presents 
a  full-PIC  model  for  the  near-thruster  plume  of  a  single  thruster.  A  multi-domain  model, 
which  splits  the  domain  into  a  near-cathode  region  and  a  near-plume  region,  is  used  to 
obtain  velocity  distribution  of  the  neutralizing  electrons.  This  simulation  is  compared  with 
one  that  assumes  a  pre-neutralized  beam  and  another  that  uses  a  floating  cathode  potential. 


I.  Introduction 

Numerical  modeling  has  been  used  extensively  to  study  the  interaction  of  electric  thruster  plume  with 
spacecraft.  Most  codes  simplify  the  plasma  dynamics  by  tracking  only  the  heavy  particles  (ions  and  neutrals), 
while  making  assumptions  about  the  distribution  of  the  electrons.  Tracking  of  electrons  using  the  full 
Particle-In-Cell  (PIC)1  algorithm  is  computationally  challenging,  due  to  a  large  difference  in  representative 
time  scales  of  the  various  species.  A  common  approach  is  to  assume  that  the  electron  distribution  responds 
to  the  ion  motion  instantaneously  according  to  the  Boltzmann  relationship,  based  on  values  of  user  supplied 
reference  parameters.  Such  a  hybrid-PIC  algorithm  cannot  correctly  resolve  the  physics  in  the  near-thruster 
region,  which  is  dominated  by  a  non-neutral  plasma  and  an  interaction  of  neutralizing  electrons  with  the 
beam.  No  simulation  models  are  currently  available  to  investigate  the  near-thruster  plume,  and  the  ion  beam 
neutralization  process  is  still  not  well  understood. 

Understanding  of  the  neutralization  process  will  become  even  more  important  for  electric  thruster  clusters 
that  are  being  considered  for  future  space  missions.  Such  a  cluster  system  may  use  a  single  cathode  to 
neutralize  ion  beams  emitted  from  multiple  thrusters.  Already,  several  cluster  configurations  were  tested 
experimentally  by  Beal  and  Hargus.2,3  However,  experimental  measurements  can  provide  only  a  limited 
amount  of  information  on  the  motion  of  the  electrons.  The  effectiveness  of  beam  neutralization  by  a  shared 
neutralizer  is  still  not  clear,  and  there  is  no  generally  accepted  optimal  configuration  for  electric  thruster 
clusters. 

This  paper  presents  results  from  a  numerical  modeling  of  neutralization  of  a  single  ion  thruster.  More 
specifically,  an  attempt  has  been  made  to  improve  results  obtained  in  a  previous  work.4  In  the  present  study, 
a  multi-domain  formulation  was  used  to  obtain  the  initial  distribution  of  electrons  near  the  cathode.  Previous 
work  used  a  floating  cathode,  since  the  primary  simulation  mesh  could  not  resolve  the  field  variations  near 
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the  cathode  tip.  Interaction  of  electrons  with  the  thruster  body  was  improved  as  well.  Electron  screening 
by  the  accelerator  grid  was  approximated  by  reflecting  electrons  that  have  collided  with  the  optics.  The 
new  model  showed  an  improvement  in  the  attained  plasma  parameters,  however  a  noticeable  discrepancy 
between  simulation  and  experimental  values  still  exists. 

II.  DRACO  ES-PIC  Code 

A  3D  plasma  simulation  code,  called  DRACO,  was  used  in  this  work.  DRACO  was  developed  within 
the  AFRL  COLISEUM  framework,  which  is  a  collection  of  modules  for  modeling  the  dynamics  of  electric 
thruster  plumes  and  their  interactions  with  spacecraft  surfaces.5 

The  DRACO  module  tracks  particles  on  a  Cartesian  mesh,  which  has  been  overlaid  by  a  secondary 
tetrahedral  mesh.6  This  secondary  mesh  allows  DRACO  to  resolve  surface  geometries  with  detail  beyond 
the  standard  ”  stair-case”  representation  typical  of  Cartesian  grids.  Surface  definition  is  specified  using  planar 
cuts  of  interface  tetrahedrons. 

The  main  COLISEUM  package  contains  support  for  loading  of  triangular  surface  meshes  from  input 
files  using  standard  formats  such  as  Ansys  or  Abaqus.  The  interface  mesh  is  generated  automatically  by 
DRACO’s  helper  module  called  VOLCAR.  The  actual  intersection  process  is  described  in  a  greater  detail  in 
Ref.  [7]. 

The  interface  cuts  are  used  to  perform  particle  surface  interactions,  and,  depending  on  the  chosen  solver, 
to  obtain  the  plasma  potential,  </>.  The  plasma  potential  is  computed  from  the  Poisson’s  equation, 

V20  =  -f  (1) 

using  the  DADI8  scheme.  In  the  above  equation,  p  is  the  charge  density  of  the  particles,  C/m3,  and  £q 
is  the  permittivity  of  free  space,  8.854  x  10-12  F/m.  The  charge  density  is  computed  from  the  individual 
contributions  of  the  ions  and  the  electrons,  p  =  q(ni  —  ne),  where  is  the  number  density  of  the  ions 
or  electrons.  In  the  PIC  method,  the  number  density  is  obtained  by  coupling  the  particles  with  the  grid 
through  particle  shape  factors , 

nk  =  E  WiS  (xi  -  xk )  (2) 

i 

where  Xk  is  the  position  of  a  grid  node,  and  is  the  specific  weight  of  the  macroparticle.  In  this  work,  the 
shape  and  size  of  the  particles  was  identical  to  the  Cartesian  cell.  This  first-order  representation  reduces 
the  simulation  noise  associated  with  the  zeroeth-order  (point  particle)  model,  while  still  allowing  a  simple 
particle-mesh  weighing  algorithm.  The  electric  field,  E,  is  then  computed  from 

V0  =  —E  (3) 

using  the  standard  centered  finite-difference  method.  Particle  velocity  is  adjusted  according  to  the  Lorentz 
force, 

dv  ->  { ->  ->\ 

m—  =  F  =  q^E-\-vx  Bj  (4) 


where  m 

=  particle  mass,  kg 

V 

=  particle  velocity,  m/s 

Q 

=  particle  charge,  C 

E 

=  electric  field,  V/m 

B 

=  magnetic  field,  T 

The  electro-static  (ES)  formulation,  implemented  by  DRACO,  assumes  that  dB/dt  =  0.  No  static  back¬ 
ground  field  was  used  in  the  current  simulation,  and  hence  the  force  acting  on  the  particles  was  simply 


The  equation  of  motion  for  the  particles  is 


=  F  =  qE 


(5) 

(6) 
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This  equation  is  integrated  numerically  along  with  eq.  5  using  the  leapfrog  method  with  a  finite  timestep  At. 
Final  position  of  the  particles  is  checked  for  surface  interactions.  Particles  leaving  the  simulation  domain 
are  either  removed  from  the  simulation,  or  are  reintroduced  according  to  prescribed  boundary  conditions. 
New  particles  are  introduced  by  sampling  particle  sources.  This  process  repeats  for  a  user  specified  time 
duration. 


III.  Neutralization  Modeling 


A.  Thruster  Model 

The  ion  thruster  used  in  this  study  is  based  on  the  40cm  NASA  Evolutionary  Xenon  Thruster  (NEXT). 
Surface  mesh  of  the  thruster  is  shown  in  figure  1.  It  should  be  noted  that  a  dimensional  drawing  of  the 
thruster  was  not  available  to  the  authors,  and  hence  the  thruster  geometry  was  generated  by  collecting  data 
from  several  relevant  sources.9,10 


(a)  x  —  y  plane 


(b)  y  —  z  plane 


Figure  1.  Surface  mesh  of  the  ion  thruster.  Large  shaded  region  indicates  source  triangles  emitting  ions. 
Electrons  are  injected  from  the  small  black  region  at  the  tip  of  the  cathode.  The  physical  curvature  of  the 
ion  optics  was  used  to  introduce  curvature  to  the  ion  beam.  The  normal  vectors  of  source  elements,  shown  by 
arrows  in  figure  (b),  roughly  indicate  the  initial  divergence  of  the  beam. 


COLISEUM  provides  a  general  support  for  surface  sources.  Various  particle  production  models  can  be 
associated  with  a  collection  of  surface  triangles.  Generally,  the  particles  are  introduced  with  a  mean  velocity 
in  the  direction  of  the  surface  normal  vector  such  that  the  physical  curvature  of  the  surface  will  result  in  a 
divergence  of  the  ion  beam.  This  feature  was  used  in  the  present  work,  as  can  be  seen  in  figure  1(b),  which 
shows  the  normal  vectors  of  the  surface  elements.  Curvature  of  the  surface  mesh  lead  to  beam  divergence  of 
approximately  15°.  Beam  flatness  (ratio  of  current  density  between  the  centerline  and  the  edge)  was  adjusted 
by  biasing  the  mass  production  rate  of  the  source  elements,  according  to  current  density  measurements  of 
Soulas.11  The  thruster  was  emitting  1.2 A  of  beam  current,  with  ions  injected  at  an  average  velocity  of 
34,400m/s  (corresponding  to  3510s  ISP)  and  a  temperature  of  O.leV.12 

Several  assumptions  were  made  about  the  plume  dynamics.  First,  collisions  were  ignored.  The  mean  free 
path,  Am,  for  an  electron- ion  collision  can  be  estimated  using13 


na 


167T£o  m2vA 
ne 4 


(7) 


For  plasma  density,  n  ~  1015  m-3  and  electron  velocities,  v  ~  106  m/s,  Am  is  0(l)m.  This  length  is  similar  to 
the  characteristic  dimension  of  the  domain.  Coulomb  collisions  thus  do  not  play  a  significant  role.  Neutral- 
neutral  and  neutral- ion  collisions  are  expected  to  be  more  frequent,  leading  to  a  particle  scatter  and  creation 
of  charge-exchange  ions.  While  the  effect  of  these  collisions  is  significant  in  the  study  of  interaction  of  the 
plume  with  the  spacecraft,  this  work  concentrated  on  the  dynamics  of  the  electrons  and  their  interaction 
with  the  beam  ions. 
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Second,  the  plume  was  assumed  to  consist  only  of  singly-charged  ions  and  electrons.  While  the  neutral 
density  near  the  thruster  exit  can  exceed  the  ion  density,  the  neutrals  interact  with  ions  exclusively  through 
collisions.  In  the  absence  of  collisions,  tracking  of  neutrals  would  only  slow  down  the  simulation.  Doubly- 
charged  ions  were  not  included  since  their  actual  distribution  was  not  known.  Furthermore,  production  of 
doubly-charged  ions  is  not  desired  because  it  leads  to  faster  thruster  erosion. 

Finally,  the  thruster  was  assumed  to  be  a  perfect  conductor.  Electrons  absorbed  by  the  thruster  shell  were 
re-injected  from  the  cathode  at  the  next  timestep.  Collection  of  an  electron  current  in  a  space  environment 
would  lead  to  a  decrease  in  the  thruster  potential,  thus  preventing  an  excessive  build-up  of  negative  charge. 
However,  the  current  version  of  DRACO  does  not  contain  functionality  to  model  this  behavior. 

The  NEXT  ion  thruster  uses  a  two-grid  ion  optics  design.  The  accelerator  grid  in  such  a  design  is  fixed  at 
potential  below  the  ambient  plasma,  specificaly  at  -257V  for  the  power  setting  used  in  this  work.11  However, 
the  magnitude  of  the  potential  drop  through  the  beamlet  centerline  is  expected  to  be  smaller,  due  to  the 
presence  of  positive  space-charge.  Numerical  modeling  of  Kafafy14  showed  potential  well  of  about  150V 
using  a  similar  configuration.  The  accelerator  grid  is  thus  responsible  for  screening  out  external  electron 
from  entering  into  the  thruster.  An  electron  would  need  to  acquire  150eV  of  kinetic  energy  to  backstream 
into  the  thruster,  which  greatly  exceeds  the  potential  difference  between  the  cathode  and  the  beam.  A 
single  finite  simulation  mesh  cannot  simultaneously  resolve  both  the  beam  and  optics  region,  due  to  large 
difference  in  associated  dimensions.  Instead,  electron  screening  was  approximated  by  elastically  reflecting 
electrons  flowing  into  the  optics. 


B.  Dimensional  Scaling  and  Boundary  Conditions 

Interaction  of  the  electrons  with  the  ions,  and  their  containment  in  the  beam  was  examined  by  modeling  a 
pre-neutralized  beam.  Beam  pre-neutralization  was  accomplished  by  injecting  both  ions  and  electrons  from 
the  optics.  Due  to  symmetry,  only  a  quarter  domain  was  simulated.  The  domain  is  shown  in  figure  2(a). 
The  cell  size  was  set  to  2  x  10-4  ^  A £>,  and  the  simulation  contained  50  x  50  x  90  cells. 


(a)  simulation  domain 


(b)  charge  density 


Figure  2.  Simulation  domain  used  to  study  electron  dynamics  in  a  pre-neutralized  beam.  Right  plot  (b)  shows 
charge  density  at  steady-state  if  a  uniform  cell  size  of  2cm  was  used.  Large  cell  size  leads  to  formation  of  a 
virtual  anode. 


The  simulation  was  performed  on  a  thruster  scaled-down  by  a  factor  of  100:1.  Plasma  density  at  the 
thruster  exit  was  retained  by  decreasing  the  emission  current  by  10,000  (100  x  100).  This  scaling  was 
necessary,  since  resolving  the  Debye  length  on  the  full-sized  domain  would  require  a  numerically  excessive 
number  of  computational  nodes  (over  1  billion).  Total  number  of  nodes  could  be  reduced  by  increasing  the 
simulation  cell  size  to  about  lOOA^.  However,  as  figure  2(b)  shows,  doing  so  results  in  a  development  of 
a  virtual  anode15  at  the  thruster  exit,  despite  the  thruster  injecting  equal  electron  and  ion  currents.  The 
virtual  anode  develops  since  the  PIC  formulation  replaces  point  sized  particles  with  particles  the  size  of  the 
cell.  No  detail  is  available  at  length  scales  smaller  than  the  cell  size.  Furthermore,  A d  specifies  the  smallest 
distance  at  which  quasi- neutrality  can  be  assumed.  Electron  motion  is  influenced  by  electric  fields  which 
arise  due  to  local  charge  non-neutralities.  A  simulation  cell  which  is  several  orders  of  magnitude  larger  than 


4  of  12 


The  29th  International  Electric  Propulsion  Conference ,  Princeton  University, 
October  31  -  November  4,  2005 

Distribution  A:  Approved  for  public  release;  distribution  unlimited. 


the  Debye  length  is  not  capable  of  resolving  these  charge  variations,  thus  proper  mixing  of  electrons  with 
the  ions  is  not  achieved. 

The  Neumann  ( d(j)/dh  =  0)  boundary  condition  was  specified  for  the  potential  solver  on  all  external 
faces.  A  reflective  boundary  condition  was  applied  for  particles  along  the  planes  of  symmetry.  Particle 
motion  through  the  remaining  domain  boundaries  was  controlled  by  an  energy  boundary  condition.  As  was 
described  in  a  greater  detail  in  Ref.  [4],  conservation  of  energy  dictates  that,  in  absence  of  other  potential 
drops  and  energy  sources,  a  particle  traveling  from  rest  through  a  potential  hump  must  reach  a  velocity 
inflection  point.  Since  kinetic  energy  must  remain  non- negative,  the  particle  motion  will  reflect,  and  the 
particle  will  be  trapped  in  a  potential  hill.  Due  to  a  limited  domain  span,  the  inflection  point  may  be  located 
outside  of  the  domain  boundary.  This  case  then  leads  to  a  removal  of  electrons  which  should  have  been 
retained  by  the  simulation.  The  energy  boundary  attempts  to  retain  these  electrons  by  reversing  velocities 
of  any  particles  with  kinetic  energy  insufficient  to  escape  the  potential  drop.  The  required  potential  energy 
is  calculated  from  the  difference  between  the  beam  core  potential  and  (j)^. 

C.  Cathode  Model 

Jameson,  et.  al.16  measured  approximately  5V  of  potential  drop  through  the  cathode,  leading  to  electron 
exit  velocities  of  approximately  106m/s.  A  cathode  with  a  keeper  exit  diameter  of  1.2cm  was  used  in  this 
work.  From  j  =  nev,  the  electron  density  at  the  cathode  exit  is  ~  6  x  1016m/s.  Electron  temperature  at  the 
exit  is  ^  leV,  yielding  ^  3  x  10-5m/s,  or  about  l/10th  of  the  cell  size.  Injection  of  electrons  from  the 
cathode  without  resolving  the  Debye  length  at  the  cathode  tip  resulted  in  development  of  a  virtual  cathode.7 
The  only  electrons  that  were  able  to  escape  were  those  born  with  a  strong  radial  velocity  component. 

In  this  work,  two  cathode  models  are  examined.  The  first  model  uses  a  floating  cathode  with  a  limiting 
value  for  charge  density  near  the  keeper  exit.  This  formulation  is  identical  to  work  done  in  Ref.  [4].  By 
floating  the  potential,  the  strong  axial  potential  gradient  is  eliminated,  and  the  electrons  can  leave  the 
cathode.  Charge  density  around  the  cathode  was  limited  to  prevent  development  of  strong  self-induced 
fields.  The  electrons  were  also  injected  with  small  velocities,  so  they  could  immediately  start  following  the 
electric  field. 


(a)  cathode  subdomain  (b)  number  density 


Figure  3.  The  small  shaded  region  indicates  the  subdomain  used  in  cathode  modeling.  Boundary  of  the  mesh 
used  in  neutralization  modeling  and  the  thruster  are  shown  for  scale.  Electrons  were  sampled  from  the  entire 
subdomain  in  subsequent  neutralization  modeling.  Plot  b)  shows  the  electron  density  and  velocity  stream  lines 
at  steady  state. 


While  this  approach  resulted  in  a  free  extraction  of  electrons,  the  final  beam  potential  and  temperature 
greatly  exceeded  experimentally  measured  values.  The  influence  of  the  cathode  model  on  the  final  results  was 
studied  by  using  a  multi-domain  formulation  to  obtain  the  initial  distribution  for  the  neutralizing  electrons. 
The  simulation  was  performed  in  two  steps,  with  first  simulation  including  only  the  near-cathode  subdomain. 
Cell  size  was  decreased  to  5  x  10-5m,  and  the  subdomain  consisted  of  30  x  40  x  60  nodes.  Potential  drop 
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of  10V  was  applied  between  the  cathode  and  the  anode,  which  was  represented  by  the  body  of  the  thruster. 
The  subdomain  is  shown  in  figure  3(a).  Boundary  of  the  primary  mesh  is  included  for  scale. 


Figure  4.  Velocity  distribution  of  the  electrons  in  the  cathode  subdomain. 

Figure  3(b)  shows  the  number  density  of  the  electrons  in  the  sub-domain  after  reaching  steady-state. 
Streamlines  of  electron  velocities  are  also  shown.  The  electrons  are  seen  to  expand  uniformly.  Geometry 
of  the  simulation  results  in  a  strong  axial  electric  field  along  the  cathode  centerline.  Electrons  born  in  this 
region  are  accelerated  out  of  the  domain.  The  remaining  electrons  are  slowed  down  by  the  potential  gradient, 
and  are  attracted  back  to  the  anode.  The  velocity  distribution  at  steady-state  is  shown  in  figure  4.  The 
electrons  retain  their  initial  injection  Maxwellian  distribution. 

The  cathode  sub-domain  served  as  a  volumetric  source  during  the  primary  neutralization  modeling. 
Position  and  velocities  of  100,000  randomly  chosen  electrons  were  sampled  to  a  data  file  at  the  end  of  the 
cathode  simulation.  The  volume  source  then  returned  a  random  entry  from  the  particle  list.  Potential  drop 
of  10V  was  retained  between  the  cathode  and  the  body  of  the  thruster,  and  charge  density  was  not  limited. 


IV.  Results 


Figure  5.  Potential,  electron  temperature,  and  charge  density  after  3  x  10  7  seconds,  for  a  single  thruster. 

Plots  in  figure  5  show  the  plasma  parameters  on  the  plane  of  symmetry  for  the  pre-neutralized  reference 
configuration  (RF).  A  distinct  beam  profile  is  seen  in  the  plot  of  the  potential.  Potential  reaches  about  4.7V 
in  the  core  near  the  thruster  exit,  and  is  seen  to  decrease  with  beam  divergence.  This  value  closely  agrees 
with  experimental  measurements.17  Electron  temperature,  computed  by  assuming  Maxwellian  distribution, 
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reaches  a  similar  value,  and  is  also  seen  to  decrease  uniformly  with  density.  Figure  5(c)  shows  the  charge 
density,  p  =  q(rii  —  ne).  A  good  neutrality  ratio  is  indicated  by  light-blue  shading.  The  charge  density  in  the 
beam  is  seen  to  be  slightly  positive,  which  leads  to  the  positive  potential  in  the  beam.  An  electron  sheath 
is  seen  to  surround  the  beam.  This  sheath  is  responsible  for  the  containment  of  the  electrons. 


(a)  electric  field,  x  component  (b)  electric  field,  z  component  (c)  electron  velocity  vectors 


Figure  6.  Electric  field  and  electron  velocity  vectors  for  the  reference  case.  Electrons  were  injected  from  the 
optics  using  Maxwellian  distribution  with  Te  =  leV. 

The  electric  field  components,  E  =  — V</>,  are  shown  in  figures  6(a)  and  6(b).  Both  the  radial  and  the 
axial  components  are  approximately  zero  in  the  bulk  of  the  beam.  Hence,  the  acceleration  of  the  electrons 
is  expected  to  be  limited  to  the  regions  near  the  edge  of  the  beam,  with  electrons  moving  at  constant 
velocities  inside  the  beam  core.  The  motion  of  the  electrons  is  highly  random  (fig.  6(c)),  even  though  they 
were  originally  injected  in  the  axial  direction,  using  a  Maxwellian  source  with  Te  =  leV.  Due  to  their  high 
mobility,  the  electrons  seem  to  have  only  a  weak  memory  of  their  injection  distributions. 


Figure  7.  Comparison  of  numerical  temperature  to  the  polytropic  model,  and  comparison  of  simulation  electron 
density  to  prediction  using  Boltzmann  model. 


Maxwellian  temperature  obtained  from  the  simulation  is  compared  to  the  polytropic  relationship 

’■-"ter 

for  three  values  of  7  in  figure  7(a).  Reference  temperature  and  density  were  chosen  to  correspond  to  the  values 
in  the  beam  core,  4.2eV  and  2.5  x  1015  m-3,  respectively.  Neither  of  the  three  chosen  gamma  values  was  able 
to  produce  an  exact  match,  however,  the  temperature  seems  to  roughly  follow  the  polytropic  relationship 
with  7  ^  1.4. 

Numerical  electron  density  was  also  compared  against  the  Boltzmann  relationship,  which  states  that  a 
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direct  relationship  exists  between  plasma  potential  and  plasma  density, 


ne  =  no  exp 


(<tl± A 

V  kT0  ) 


(9) 


Again,  plasma  properties  in  the  beam  core  were  used  for  the  reference  values.  Reference  potential  was  set 
to  4.7V.  The  relationship  was  computed  using  both  constant  reference  temperature  (4.2eV),  and  polytropic 
temperature  with  7  =  1.4.  Generally,  the  agreement  is  not  very  good,  as  figure  7(b)  shows.  Best  agreement  is 
achieved  near  the  core,  which  is  expected,  since  this  location  corresponds  to  the  point  at  which  the  reference 
values  were  sampled.  The  simulation  electron  density  drops  off  faster  than  predicted  by  the  model.  The 
disagreement  is  reduced  by  using  the  polytropic  temperature,  however,  a  significant  discrepancy  still  remains. 
Adjustment  of  the  reference  parameters  would  result  in  a  better  agreement  of  the  Boltzmann  model  with 
the  numerical  electron  distribution;  this  approach  however  requires  prior  knowledge  of  the  solution. 


B.  Single  Thruster 

1.  Floating  Cathode 


The  simulation  domain  used  in  the  single  thruster  neutralization  modeling  can  be  seen  in  figure  3(a).  Due 
to  symmetry,  only  a  half  domain  was  simulated.  Reflective  particle  boundary  condition  was  used  along 
the  symmetric  face.  The  grid  dimensions  were  50  x  100  x  90,  and  the  mesh  used  a  uniform  cell  spacing  of 
2  x  10_4m. 


(a)  potential 


(b)  electron  temperature 


(c)  charge  density 


Figure  8.  Potential,  electron  temperature,  and  charge  density  after  3  x  10  7  seconds,  for  a  single  thruster, 
floating  cathode. 


Figure  strip  8  shows  potential,  temperature  and  charge  density  obtained  by  floating  the  potential  on  the 
cathode.  Results  for  this  configuration  (NSF)  were  expected  to  agree  with  the  reference  case  RF,  but  a  quick 
comparison  with  figure  8  indicates  a  high  degree  of  divergence.  The  beam  shape  is  no  longer  well  resolved. 
Furthermore,  the  beam  potential  has  increased  to  27V  and  the  maximum  electron  temperature  has  increased 
to  34eV. 

The  high  values  of  beam  potential  and  temperature  point  to  a  non-neutral  beam.  However,  over¬ 
saturation  of  the  beam  with  electrons  by  increasing  the  cathode  current  led  to  turning  of  the  beam  towards 
the  cathode,  but  the  potential  drop  across  the  beam  did  not  change  significantly.  The  beam  temperature 
also  remained  high. 

2.  Multi-domain  Cathode  model 

Therefore,  lack  of  electrons  did  not  seem  to  be  the  main  factor  contributing  to  the  non-physical  results.  The 
influence  of  the  cathode  model  on  the  results  was  investigated  by  replacing  the  floating  potential  model  with 
the  multi-domain  formulation.  Figure  9  shows  plasma  properties  obtained  with  the  new  model.  Important 
to  note  is  the  region  near  the  cathode.  The  charge  density  plot  shows  a  clear  turning  of  the  electron 
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(a)  potential 


T  .6-  (eV) 


(b)  electron  temperature 


(c)  charge  density 


Figure  9.  Potential,  electron  temperature,  and  charge  density  after  3  X  10  7  seconds,  for  a  single  thruster. 


plume  towards  the  beam.  This  behavior  was  not  achieved  in  the  floating  cathode  case.  The  electron  sheath 
surrounding  the  beam  is  also  better  resolved,  and  the  overall  range  of  p  has  decreased  which  indicates  a 
better  neutrality  ratio.  Beam  potential  remains  above  the  measured  values,  however,  the  peak  magnitude 
has  decreased  to  about  23V.  The  beam  profile  is  also  free  of  any  significant  asymetrical  anomolies.  Similarly, 
the  electron  temperature  remains  high,  but  the  temperature  decay  became  more  polytropic. 

C.  Analysis 


nd.e-  (#/m3) 

I1.3E+18 
2.0E+15 
1.5E+15 
1  .OE+15 
5.0E+14 
2.4E+14 
0.0E+00 


nd.e-  (#/m3) 

S1.4E+15 
1. OE+15 
5.0E+14 
2.4E+14 
0.0E+00 


(a)  RF 


(b)  NSF 


(c)  NSM 


Figure  10.  Electron  density  contour  for  the  three  cases.  A  clear  distinction  is  seen  between  the  reference  case 
and  the  two  cases  in  which  electrons  were  injected  from  the  cathode,  multi-domain  cathode  model. 

Although  the  new  model  lead  to  some  improvements  in  the  results,  it  seems  that  it  failed  to  resolve  a 
fundamental  problem  existing  in  the  simulation.  Figure  10  shows  the  electron  densities  for  the  three  cases. 
A  clear  difference  is  seen  between  the  reference  case,  and  the  two  cathode  cases.  The  electron  density  in  the 
reference  case  follows  the  Boltzmann  relationship,  with  highest  density  coinciding  with  the  beam  core.  The 
electrons  instead  seem  to  concentrate  along  the  edges  of  the  beam  in  the  two  cathode  runs. 

Existence  of  an  almost  uniform  high-temperature  region  along  beam  axis  in  the  NSF  case  indicates  mixing 
of  electron  streams  with  opposing  directions.  Thus,  the  collective  dynamics  of  the  electrons  seem  to  be  driven 
by  oscillations  around  the  beam  core.  In  other  words,  using  a  1-D  approximation,  electrons  injected  at  the 
cathode  fall  into  the  potential  well  created  by  the  beam.  The  velocity  of  the  electrons  increases  until  they 
pass  the  beam  centerline.  The  velocity  then  begins  to  decrease  as  they  travel  up  the  well.  The  electrons 
come  to  a  stop  at  a  point  where  all  kinetic  energy  has  been  exhausted.  Assuming  initial  injection  at  Om/s, 
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the  potential  drop  on  both  sides  of  the  well  will  be  equal.  The  increased  electron  densities  seen  along  the 
beam  edge  in  cases  NSF  and  NSM  are  a  direct  consequence  of  electrons  coming  to  a  stop  and  turning  around. 


Figure  11.  Electron  and  ion  velocity  distribution  at  the  end  of  simulation. 

The  potential  in  the  beam  adjusts  as  a  direct  response  to  the  oscillations  of  the  electrons.  Retention  of 
the  electrons  in  the  beam  requires  that  the  PEbearn  >  KEe:Tnax.  Using  PE  =  KE  and  total  potential  drops 
between  the  cathode  and  the  beam  of  4.7V  and  25V,  the  mean  velocity  of  the  electrons  oscillating  around 
the  beam  core  is  expected  to  be  1.3  x  106  and  2.9  x  106  m/s  for  cases  RF  and  NSx,  respectively.  This  is 
confirmed  by  the  tail  of  electron  velocity  histograms  shown  in  figure  11(a).  The  dotted  line  (NSF)  shows  the 
distribution  for  the  original  floating  cathode  model.  Comparison  with  the  reference  curve  shows  an  increase 
in  mean  velocity.  Even  more  important  is  the  widening  of  the  curve,  which  is  indicative  of  a  temperature 
increase.  The  shape  of  the  curve  remains  close  to  Maxwellian.  The  dashed  line  (NSM)  indicates  the  electron 
velocity  obtained  with  the  new  model.  A  small  increase  in  the  number  of  fast  moving  particles  is  seen, 
but  overall,  the  distribution  remains  comparable  to  NSF.  A  greater  difference  is  seen  in  the  ion  velocities. 
Decrease  in  the  size  of  the  high  potential  region  in  the  NSM  case  results  in  fewer  electrons  being  slowed 
down  by  the  strong  potential  gradient.  The  mean  beam  velocity  for  case  NSM  is  closer  to  the  reference 
configuration. 

In  order  to  obtain  the  Boltzmann  relationship,  the  amplitude  of  the  electron  oscillations  must  decrease 
with  time.  This  is  analogous  to  the  classical  example  of  stability,  in  which  a  ball  has  been  placed  into  a 
spherical  cup.  Displacement  of  the  ball  from  the  rest  position  at  the  bottom  will  result  in  simple  oscillations 
about  the  bottom  of  the  cup.  In  the  absence  of  dissipative  forces,  the  amplitude  of  oscillations  will  not 
change  with  time.  However,  in  a  realistic  configuration,  dissipation  of  energy  due  to  non- conservative  forces 
will  result  in  a  gradual  decrease  of  the  amplitude.  After  some  time,  the  ball  will  come  back  to  rest. 

The  simulation  presented  here  does  not  contain  any  such  dissipative  force.  Instead,  the  electrons  keep 
oscillating  around  the  potential  drop  of  the  beam.  Hence,  the  electron  density  is  not  able  to  reach  a 
Boltzmann-like  relationship.  Electrons  in  the  reference  case  are  not  strongly  influenced  by  this  simplification, 
since  they  are  born  at  the  bottom  of  the  potential  well.  Furthermore,  their  initial  velocity  is  coincident  with 
the  direction  of  the  beam  ions.  However,  the  electrons  born  at  the  cathode  originate  at  the  top  of  the 
potential  well  and  flow  into  the  beam  radially. 

Exchange  of  kinetic  energy  with  massive  particles  would  result  in  a  large  decrease  in  electron  velocities, 
while  only  slightly  influencing  the  motion  of  the  ions.  However,  as  was  mentioned  previously,  collisions  were 
ignored  in  this  case,  due  to  large  mean-free  paths,  and  hence  low  collision  frequencies.  More  probable  is  the 
transfer  of  kinetic  energy  from  the  electrons  to  plasma  waves.  Figure  12  shows  the  total  field  energy  versus 
simulation  time.  Clear  oscillations  develop  after  t  =  1.5  x  10-7  seconds.  The  meaning  of  these  oscillations 
is  still  not  clear,  however  a  decay  is  seen  around  t  =  2.4  x  10-7.  It  is  possible  that  the  numerical  setup  of 
the  problem  is  preventing  development  or  growth  of  energy  dissipating  waves. 
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Figure  12.  Total  simulation  field  energy  versus  simulation  time.  A  high-frequency  wave  is  seen  to  develop 
after  t  =  1.5  X  10— 7  seconds. 

V.  Conclusion 

A  new  simulation  model  to  study  ion  beam  neutralization  is  being  developed.  This  model  uses  a  fully- 
kinetic  formulation  in  which  both  electrons  and  ions  are  tracked  as  macro-particles.  This  formulation  avoids 
problems  associated  with  fluid  modeling  of  the  electrons,  but  introduces  numerical  difficulties.  Most  im¬ 
portant  is  the  necessity  to  resolve  the  local  Debye  length,  otherwise  the  electrons  fail  to  mix  with  the  ions. 
Computationally  excessive  number  of  nodes  would  be  required  to  resolve  the  Debye  length  on  a  full-scale 
geometry.  Instead,  a  dimensional  scaling  approach  was  used,  with  thruster  current  adjusted  such  that  plasma 
density  at  the  thruster  exit  remained  unchanged.  This  approach  allowed  for  examining  the  neutralization 
process.  Outflux  of  electrons  at  boundaries  was  prevented  by  reflecting  all  electrons  with  kinetic  energies 
insufficient  to  escape  the  potential  drop  of  the  beam. 

This  simulation  approach  was  used  to  model  neutralization  of  the  NASA  NEXT  ions  thruster.  A  reference 
case  was  setup  by  injecting  both  electrons  and  ions  from  the  optics.  The  potential  solution  showed  a  clear 
beam  profile,  with  maximum  potential  of  4.7V.  The  electron  temperature  reached  about  5eV  in  the  core, 
and  decreased  polytropically  with  density.  These  results  agree  well  with  experimental  data.  The  electron 
density  was  also  compared  to  the  Boltzmann  model,  but  the  two  curves  diverged  for  the  chosen  coefficients. 

Simulation  of  a  single  thruster  neutralized  by  an  external  cathode  was  also  studied.  The  potential  around 
the  cathode  could  not  be  resolved  correctly  using  the  primary  mesh,  due  to  the  cathode’s  small  size  and 
a  high  electron  density  near  the  tip.  Instead,  two  approaches  were  investigated.  In  the  first  model,  the 
cathode  potential  was  allowed  to  float  and  charge  density  at  the  tip  was  fixed.  The  second  model  used  a 
multi-domain  formulation.  Simulation  of  the  near  cathode  region  was  performed  first  using  a  fine  mesh, 
followed  by  sampling  of  random  electrons  to  a  data  file.  This  distribution  list  was  then  used  by  primary 
simulation  to  introduce  electrons  from  the  volume  described  by  the  fine  mesh. 

Plasma  properties  in  either  cathode  run  did  not  agree  with  the  reference  case.  The  floating  cathode 
model  resulted  in  beam  potential  of  about  27V  and  temperature  distribution  which  no  longer  followed  the 
polytropic  relationship.  The  new  cathode  model  showed  an  improvement  in  the  results,  with  beam  potential 
decreasing  to  about  23V  and  temperature  assuming  a  more  polytropic  decay.  However,  the  discrepancy 
between  results  and  experimental  measurements  remained  significant.  Closer  inspection  of  the  simulation 
results  indicates  that  introduction  of  electrons  from  the  cathode  results  in  oscillations  around  the  beam  core. 
Attaining  a  Boltzmann-like  density  requires  the  oscillations  to  decay  with  time,  however,  this  does  not  seem 
to  be  the  case  in  the  current  model.  An  investigation  of  the  primary  decay  mechanism,  and  its  numerical 
implementation,  remains  as  future  work. 
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